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ABSTRACT 

Classical Sweet-Parker models of reconnection predict that reconnection rates depend inversely on 
the resistivity, usually parameterized using the dimensionless Lundquist number (S). We describe 
magnetohydrodynamic (MHD) simulations using a static, nested grid that show the development of 
a three-dimensional instability in the plane of a current sheet between reversing field lines without a 
guide field. The instability leads to rapid reconnection of magnetic field lines at a rate independent 
of S over at least the range 3.2 x 10 3 < S < 3.2 x 10 5 resolved by the simulations. We find that 
this instability occurs even for cases with S < 10 4 that in our models appear stable to the recently 
described, two-dimensional, plasmoid instability. Our results suggest that three-dimensional, MHD 
processes alone produce fast (resistivity independent) reconnection without recourse to kinetic effects 
or external turbulence. The unstable reconnection layers provide a self-consistent environment in 
which the extensively studied turbulent reconnection process can occur. 


1. INTRODUCTION 

During magnetic reconnection, magnetic field lines 
change topology, resulting in the conversion of magnetic 
energy into both thermal energy and kinetic energy of 
bulk flows and non-thermal particles. The rate at which 
this process occurs in the classical Sweet-Parker picture 
(Sweet |[l958 |Parker|[l957 ) depends on the Lundquist or 


magnetic Reynolds number S = vaL/tj , where va is the 
Alfven speed, L a characteristic length of the system, and 
ri the resistivity. The Sweet-Parker rate is orders of mag¬ 
nitude too slow to explain the fast reconnection seen in 
low resistivity plasm as during solar flares and sawtooth 
crashes in tokamaks (Yamada et al.||2010). Because it is 


a fundamental plasma process, reconnection is thought 
to be important in astr ophysical environments a s diverse 
as the helio sphere (e.g. Edmon dson et al.||2010 1 and mi¬ 


croquasars (Khia lTet al.|20l5 ). 


l(Biskam 

p 1986 

Loureiro et al. 

2007 

Huang & Bhat- 

tacharjeeJ2(J13f, a super-Alivenic, smai 
has provided a mechanism to greatly 

[-scale instability, 
speed up Sweet- 


Parker reconnection. However, this instability has pri¬ 
marily been studied in two dimensions (2D), assuming 
symmetry in the plane of the current sheet. The reason 
for this dimensional reduction is that the reconnection 
process is inherently multi-scale, with a large separation 
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between the global scale of the reconnection layer and 
the resistive length where the instability grows. Even 2D 
simulations tax state of the art computational resources 
if unif orm grids are use d. 


Lazarian & Vishniac (1999) argued that reconnection 
in the presence of any sort of turbulence would be fast, 
because the turbulence would drive many points of con¬ 
tact between the opposed field lines. This idea has been 
put on a rigorous mathematical basis ( Eymlc_eta^20lO 


as rev iewed by Lazarian et al. (2015a I and Lazarian et al. 
(|2015||. Indeed, recent modelling suggests that turbu- 
lent reconnection may be responsible for the rad io and 


gamm a-ray emission from accreting black holes (Singh 


et al. 2015). However, these ideas all require that re¬ 


connection proceed at a large fraction of va- Numerical 
models examining reconnection in f orced turbul ence sup¬ 
port this theory, starting with (Kowal et al. 2009). In 
this work, we demonstrate that turbulent reconnection 
proceeds in a very similar fashion when the turbulence 
is self-generated from an instability of the reconnection 
layer itself. 

Here, we describe a set of nested grid simulations that 
model the reconnection layer in 3D over a broad range 
of S, without any forcing or guide field. These simula¬ 
tions show that a startlingly fast, 3D, instability occurs 
in the plane of the current sheet, which was assumed 
uniform in the 2D simulations. This instability drives a 
large increase in the rate of reconnection, that we show 
remains independent of S over two orders of magnitude 
of variation in the resistivity. 
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Boozer (2012b 2013) has argued that reconnection re¬ 


quires a point geometry to proceed, so that it is an in¬ 
herently three-dimensional (3D) process. The work we 
describe here demonstrates that such a 3D geometry nat¬ 
urally arises even from 2D initial conditions, resulting in 
fast reconnection apparently independent of S. 

Previous work in this held has shown 3D instability, 
but has not provided a clear dem onstration of indepen¬ 


dence of reconnection rate from S. Dahlburg et al. (2003 
20051 focused on the case of a current sheet with a strong 


guide held, and found a 3D instability set in for a weak 
enough guide held, which they called a secondary insta¬ 
bility. However, they did not measure the sc aling of the 
reconnection rate with S. Lapenta & Bettarini (2011) re¬ 
ported the breakdown of an initially 2D Harris sheet into 
a fully 3D reconnection region with greatly enhanced re¬ 
connection rate. An MHD kink instability on a central 
plasmoid was followed by a Rayleigh-Taylor instability 
driven by the reconnection jet interacting with the plas- 
moids at the ends of the layer. However, again, no test 
of the dependence on S was performed. Another numer¬ 
ical experiment has shown that thin, 3D, current layers 
are unstable to infinitesimal perturbations and reconnect 
at a rate ap p aren tly independent of Lundquist number 
S (Beresnyak||2013 1, but only a factor of three variation 


in S was explored. Edmondson et al. 


(2010) studied the 


formation of coronal current sheets due to photospheric 
forcing in a global, 3D, AMR simulation. They concluded 
that the dynamics of the current sheet were 3D, allowing 
a steady rather than the bursty reconnection rate found 
by 2D models of the plasmoid instability. Finally, 3D re- 
connection in the collis ionle ss limit has been explored by 
Daughton et al. (2011) and Pritchett (2013|). That work 
largely focused on particular kinetic effects that drive 
dissipation at the smallest scales. 

In Sect. [2] we describe our computational methods, 
while in Sect. [3] and [4] we present our results. Finally, 
we discuss the implications in Sect. [5| 

2. METHODS 

We use the mesh refinement code Enzo , which solves 
t he com pressible, adiabatic, resistive, MHD equations 
(Bryan et al. 20141. We use static refinement to focus 
computationareffort on the current sheets. Our compu¬ 
tations are performed within a cubic, 3D volume, with 
periodic boundary conditions in all three dimensions. 
From the available algorithmic options, we choose piece- 
wise linear reconstruction, the HLLD Riemann shock¬ 


capturing so lver, and constrained transport (Gardiner & 
Stone 2005) t o ensu re V • B = 0 to machine precision 
(ICollins et al.| |2010[). We perfo rmed all analysis using 
the yt toolkit ( Turk et al.||2011 ). 

Our initial condition is a pair of oppositely directed, 
parallel current sheets to accommodate the periodic 
boundary conditions, each perturbed foll owin g the GEM 
Reconnection Challenge (Birn et al. 2001) to initiate 
Sweet-Parker reconnection. All but two oi our runs are 
initialized with low-amplitude, 3D velocity perturbations 
with mean Alfven Mach number ( v/va ) ~ 4.3 x 10~ 5 . 
These perturbations have a spectrum Vk oc k ~ 4 with 
wavenumbers ranging from k m i n < k < k max . We choose 
kmin/ 27t = 10 and k max /2n — 15, except for run C+, 
which has k m i n /2TT = 30 and k max /2Ti = 35. We do not 
continue to force the velocity field during the simulation. 


We normalize all lengths to the size of the box L = 1, 
densities to the initial density at the center of the sheets 
po = 1, and times to the Alfven crossing time of each 
sheet tA = Sq/va where va = Ho/V^Too — 3.2 is the 
Alfven speed of the upstream plasma. S 0 = 0.02 is the 
scale length of the initial current sheet. Table [l] lists the 
parameters of our runs. 

A minimum resolution requirement for reconnection is 
the proper resolution of the Sweet-Parker current layer, 
whose width 5$p — Lj\/~S , where L is the length of the 
layeiQ We define two grid refinement regions covering 
the entire plane of the current sheet (0 < [ x r , y r ] < 1) 
with a height z r ~ 12.5^o centered on each of the initial 
current sheets. These refined regions have two levels of 
refinement atop a 128 3 base grid, leading to an effective 
resolution of 512 3 in the current sheet centers, except 
for run A*, which has three levels of refinement, for an 
effective 1024 3 resolution. 

Figure [T] shows Sgp of the initial Sweet-Parker current 
sheet at t = 75tA, long before any unstable perturbations 
have grown to significant amplitudes. All runs with S < 
10 5 have current sheet widths that agree well with the 
Sweet-Parker prediction, because they are resolved by 
10 zones across the sheets. The run with S = 3.2 x 10 5 
demonstrates the effects of marginal resolution, while the 
run with S = 3.2 x 10 6 is only resolved by ~ 3 zones, 
and is a factor of four too thick. We do not use this last 
run (run J) in our subsequent analysis, although it serves 
as an important limit on the numerical resistivity of our 
code. 


Table 1 

Run data 


run 

s 1 



v 2 

7 3 


notes 

A 

3.2 

X 

10 5 

10“ 5 

-3.2 x 

10" 3 


A* 

3.2 

X 

10 5 

10“ 5 

-3.3 x 

10~ 3 

double resolution 

B 

3.2 

X 

10 4 

10 -4 

-5.6 x 

10" 3 


C 

1.6 

X 

10 4 

2 x 10 -4 

-4.8 x 

10" 3 


C+ 

1.6 

X 

10 4 

2 x 10 -4 

-1.8 x 

10" 3 

kmin/ 27r = 30 perturbation 

D 

8.0 

X 

10 3 

4 x 10 -4 

-2.1 x 

10“ 3 


E 

3.2 

X 

10 3 

10 -3 

-1.4 x 

10“ 3 


F 

3.2 

X 

10 2 

10 -2 

- 


stable to 3D instability 

G 

3.2 

X 

10 3 

10 -3 

- 


no initial perturbations 

H 

3.2 

X 

10 5 

10 -5 

- 


no initial perturbations 

J 

3.2 

X 

10 6 

10-® 

- 


underresolved, unanalysed 

Lundquist 

number 2 Resistivity in 

code 

units 3 Decay rate of 


magnetic energy (see text) 


3. FIELD DYNAMICS 

Reconnection in our models begins at the Sweet-Parker 
rate expected for a stable field reversal, as shown by the 
width of the current sheet. This leads to the initial slow 
decline of the volume integrated magnetic energy for all 
simulations (Figure [2]), as well as the low values of in¬ 
tegrated kinetic energy. Instability along the plane of 
the current sheet then sets in, driving far faster recon¬ 
nection, and transferring energy from the magnetic field 
into the flow, as shown by the sudden drop in magnetic 
energy and the corresponding rise in kinetic energy. The 
morphology of the onset and growth of the instability is 

1 A popular alternative is to use 5 to represent the half-width 
of the current sheet, but in that case, L is the half-length as well. 
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Figure 1. Width Ssp in units of box size of the current sheet 
during quiescent reconnection at t = 75 1^, prior to the onset of 
instability. The circles show simulations with varying Lundquist 
number S, the solid line gives the Sweet-Parker scaling, and the 
triangle shows run A* at double resolution. The dot-dashed line 
shows a resolution of 10 zones for standard (512 3 -equivalent) runs, 
while the dotted line shows 10 zones for our high resolution (1024 3 - 
equivalent) run. 


shown in the middle panels of Figure [3] while the final 
panel shows its saturated state. 

The development of the 3D instability results in the 
buckling of the current sheet in the y—z plane, with a 
characteristic wavenumber k z /2i r ~ 12 ( third pa nel of 
Fig. |3|. The simulations of Lapenta & Bettarini (2011) 
can be seen to show similar behavior, though it is not 
emphasized in their paper. They used a thin box in the 
third dimension, so they only had two wavelengths in 
that direction. Ours is a factor of six deeper in the z di¬ 
rection than theirtj^J Thus, our results give a wavenum¬ 
ber consistent with that shown in their figures. 

We find that resistivity stabilizes the 3D instability for 
S < 10 3 . Run D with S = 3.2 x 10 3 shows the instability 
clearly through the growth of kinetic energy, although 
the total amount of reconnection driven by the turbu¬ 
lence is small, because the large resistivity has already 
allowed significant laminar reconnection to occur. 

When S > 10 4 , the transition from laminar to turbu¬ 
lent reconnection begins with the rapid growth of kinetic 
energy apparently driven by a kink-type instability along 
the plasmoids in the ^-direction. However, this is not a 
classical kink instability, as the interior of the plasmoids 
is a demagnetized reconnection region, rather than a col¬ 
umn of current-carrying plasma surrounded by vacuum 
as would be true in the classical case. The growth of ki¬ 
netic energy and decay of magnetic energy must be due 
to reconnection occurring where field lines are driven to¬ 
gether by these instabilities rather than a simple rear¬ 
rangement of the horizontal field by them. 

It appears from our models that the 3D instability 


may actually grow independently of the 
stability. Run E with S < 10 4 lacks 
growth of the 2D plasmoid instability 
Lundquist numbers by ourselves and 

2D plasmoid in¬ 
evidence for the 
seen at higher 
previous authors 

(Loureiro et al. 2007; Samtaney et al. 

2009 Uzdensky 

et al. 2010: Huang & Bhattacharjee 2010: Loureiro et al. 

2012 Huang & Bhattacharjee 2013), 

but nonetheless 
albeit with a de- 

the 2D plasmoid 
r-Meshkov insta- 

) box is oriented so 

shows the growth of the 3D instability, 
lay in onset of rapid growth (Fig. p]). 

This delay occurs because growth of 
instability triggers secondary Richtmye 

2 Note that the Lapenta & Bettarini] (|2011 
their y axis corresponds to our z axis 






Figure 2. (upper) Volume integrated kinetic and magnetic en¬ 
ergies in the simulation domain as a function of time for several 
values of Lundquist number S. Letters in the legend give the run 
name in Table [l] and the gray tickmarks on the central axis give 
the times of the four panels in Figure [3] (lower) Current sheet 
width A as a function of time for run A*. The light line shows 
a linear fit to the unstable period, giving a layer growth rate of 
dA/dt ~ 3 X 10 _3 ua- 


bility (Richtmyer 1960| Meshkov 1969) that accelerates 
onset of the 3D instability, but is not required for 3D in¬ 
stability to occur. The acceleration of a flow across the 
density contrast between the plasmoid and the surround¬ 
ing current sheet along the x axis drives the Richtmyer- 
Meshkov instability (analogous to the Rayleigh-Taylor 
instability that occurs in a gravitational field). This in¬ 
stability begins when weak transient shocks from the 
formation of the plasmoids at the sides of the domain 
(x = —0.5, 0.5) pass over the density gradient produced 
by the first plasmoid to form around the initial X line at 
(x,y) = (0,0.5). The secondary instability drives weak, 
initial mixing along the X line. 

We examined three numerical issues with further runs. 
First, to determine if resolution affects our major result, 
we run our marginally resolved run A at twice the reso¬ 
lution (run A*). This run results in essentially identical 
growth rate 7 . 

Second, we checked that the instability is not purely 
numerical by performing Runs G & H, identical to the 
unstable Runs E & B respectively, except without any 
initial perturbations. The instability did not grow in ei¬ 
ther of these runs. In Run G, which has S = 3.2 x 10 3 
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and is thus stable to the 2D plasmoid instability, the 
entire 3D volume settled down into a steady, laminar re¬ 
connection at the Sweet-Parker rate, with a sheet width 
S ~ Ssp = 1-72 x 10“ 2 L. Run H, on the other hand, 
is unperturbed in the y-z plane, but is unstable to the 
plasmoid instability. In this case, we find a vigorous plas¬ 
moid instability that is entirely symmetric along the 2 
axis, demonstrating that our code is sufficiently stable 
to recover the 2D results in 3D if there are no explicit 
3D perturbations. 

Third, run C+ was performed at standard resolution 
with smaller scale perturbations (fc min /2-7T = 30). We 
find the growth rate is higher by a factor ~ 3 for the 
lower k perturbations, suggesting a wavenumber depen¬ 
dence for the underlying instability. We will pursue a 
formal stability analysis of the instability in a separate 
paper, and this wavenumber dependence represents an 
important test for that work. 


4. SCALING 

Figures [2] and [3] show three phases of reconnection: the 
slow, Sweet-Parker phase while linear instabilities grow, 
a rapid exponential phase in which 3D effects dominate 
reconnection, and finally a saturated, MHD turbulent 
phase. The Sweet-Parker reconnection rate is oc S -1 / 2 , 
implying far slower than observed reconnection at the 
large values of S typical of Solar and space plasmas. Fig¬ 
ure [4] shows the decay rate 7 during the rapid reconnec¬ 
tion phase as a function of S. While S runs over two 
orders of magnitude, 7 varies by a factor of only about 
three, with no discernible functional relationship to S. 
Thus, the 3D instability offers a fast reconnection mech¬ 
anism that occurs at a rate apparently independent of 
S, without appeal to either kinetic effects or anomalous 
resistivity. 

Once the instability saturates, the current sheet thick¬ 
ens considerably (see the right panel of Fig. [ 3 ]), and the 
picture of a steady flow of fresh field from upstream (i.e. 
from the y direction above and below the layer) no longer 
holds in our simulations. At the end of our simulation, 
there is still plenty of field left to reconnect. The tur¬ 
bulence self-consistently driven in the r eco n nection l ayer 
allows stochastic reconnection to o ccur (Lazarian & Vish- 
niac|1999 Eyinketal. 2011 2013). The thickening of the 
reconnection layer we see is consistent with their model: 
the diffusion of the large-scale magnetic fields controls 
the ultimate reconnection rate. Our per iodic boundary 
conditions, similar to those of Beresnyak (2013), do not 
allow for a self-consistent steady state to occur. How¬ 
ever, during the rapid, resistivity-independent phase of 
the instability, the reconnection region grows. 


5. DISCUSSION 

We have shown that for Lundquist numbers S > 3.2 x 
10 3 , current sheets become turbulent in the direction per¬ 
pendicular to the field along the sheet, leading to fast, 
3D, magnetic reconnection (i.e. decay rate of magnetic 
energy independent of resistivity y). At high S, indi¬ 
vidual 2D plasmoids rapidly lose their identities, as the 
current sheet splits into filaments parallel to the field di¬ 
rection. These 3D effects, driven by rapidly growing in¬ 
stabilities along current sheets, appear essential to under- 


standing reconnection. This provides strong sui 

oport to 

the geometric ideas advanced by 

Boozer 

(2012b 

a, 

2013) 


as well as the turbulent reconnectio n model developed b y 
Lazar ian, Vishniac, a nd coworkers (jLazarian fe Vishniac 


1999[ Lazarian et al. 2015 a). To demonstrate the latter 


point, in the bottom panel of Figure[2]we show the recon¬ 
nection layer width A as a function of time. During the 
rapid growth phase, A oc t, in pleasing agreem ent with 
the theory described in Lazarian et al. (2015a). We find 
that the layer expansion velocity is dA/dt ~3x lO -3 ^, 
roughly a f actor of 4-5 smaller than that reported by 


Beresnyak (2013). We suspect the discrepancy is due 
to theTact that his simulations include a guide field, al¬ 
though the lower diffusivity of his pseudo-spectral code 
could also play a role. Nevertheless, the agreement on the 
linear form of the growth rate despite our different setups 
and codes supports the turbulent reconnection model. 

As a result of the 3D instability, the initial current 
sheet develops into a thick region of MHD turbulence. 
Lazarian et al. (|2015al speculate that the growth in the 


plane perpendicular to the reconnection (i.e. along the 
z direction in our simulations) could be due to Kelvin- 
Helmholz instability. Once the turbulent state is reached, 
the decay of magnetic energy in our model slows dramat¬ 
ically, back to a rate comparable to the initial Sweet- 
Parker rate. However, we stress that the state of the 
system after the rapidly reconnecting, unstable phase is 
radically different from the state before it, and the slow 
subsequent reconnection may depend on the geometry 
we chose for our simulations, which does not continue to 
force t he system on l arge scales, unlike, for example Solar 
flares Dudfk et al. (2014). However, this turbulence nat¬ 
urally proHuces - tn£)conditions_£eguiredJorfast_stochas^ 


tic reconnection (Lazarian fe Vishniac|1999; Eyink et al. 


2011 2013). 

Future work employing deeper grid hierarchies, adap¬ 
tive resolution elements placed in regions of high J, 
and simulations with large scale forcing through inflow 
boundary conditions will clarify the outcomes of the in¬ 
stability and its role in understanding reconnection in as- 
trophysical environments. Specifically, by removing pe¬ 
riodic boundary conditions (or isolating them via deep 
AMR hierarchies), we will be able to m ake a more de¬ 
tailed stud y of a key prediction of the |Lazarian Sz VishT] 
niac (1999) model, the rate of broadening of reconnection 
layer. 
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Figure 3. Slices of current density |J| at four times in run A, with S = 3.2 x 10 5 . Each panel shows the lower current sheet in our 
simulation box (the upper sheet looks morphologically similar at each time). Note that the instability grows in both the y—z and x—z planes 
at roughly the same time. 
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